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The eigenmodes and the vibrational density of states of the ground state configuration of graphene 
clusters are calculated using atomistic simulations. The modified Brenner potential is used to 
describe the carbon-carbon interaction and carbon-hydrogen interaction in case of H-passivated 
edges. For a given configuration of the C-atoms the eigenvectors and eigenfrequencies of the normal 
modes are obtained after diagonalisation of the dynamical matrix whose elements are the second 
derivative of the potential energy. The compressional and shear properties are obtained from the 
divergence and rotation of the velocity field. For symmetric and defective clusters with pentagon 
arrangement on the edge, the highest frequency modes are shear modes. The specific heat of the 
clusters is also calculated within the harmonic approximation and the convergence to the result for 
bulk graphene is investigated. 

PACS numbers: 63.22.Kn, 63.22.Rc, 65.80.Ck 



I. INTRODUCTION 

During the past few years considerable attention has 
been paid to the study of the spectral properties such 
as the energy spectrum, the eigenmodes (i.e. phonon 
spectrum), and the phonon density of states of car- 
bon allotropes such as graphene^ ^, carbon nanotubeP^, 
fullereneJ^M^, graphene flakeJ^ and carbon cluster J^^. 

Raman spectroscopy is a widely used tool to study the 
vibrational modes of carbon based materials^^ Be- 
cause the normal modes of molecules are unique, they 
have their own set of Raman frequencies. Similarly for 
graphene clusters their phonon modes will give us infor- 
mation on e.g. their structural and dynamical properties. 

Lattice dynamics calculations of a single graphene 
sheet were done by Rao et based on C-C force 

constanti^^^^ which were obtained from a fit to the 
two dimensional, experimental phonon dispersion. Good 
agreement with Raman spectra of single wall car- 
bon nanotubes was found. Tommasini et al.^^ pre- 
sented a semiempirical approach for modeling the off- 
resonance Raman scattering of molecular models of con- 
fined graphene and compared with the results from den- 
sity functional theory calculations. Pimenta et al.^^ re- 
viewed the defect-induced Raman spectra of graphitic 
materials and presented Raman results on nanographites 
and graphenes. Ryabenko et al.^^ analyzed the spectral 
properties of single-wall carbon nanotubes encapsulat- 
ing fullerenes and found that fullerene alters the aver- 
age diameter of the electron cloud around the single wall 
carbon nanotube. Michel et al.^ presented a unified 
theory of the phonon dispersion and elastic properties 
of graphene, graphite and graphene multilayer systems. 
They started from a fifth-nearest-neighbor force-constant 
model derived from the full in-plane phonon dispersions 
of graphitd^, and extended this model with interplanar 
interactions, and found that the graphite eigenfrequency 
^B2g^ ^127 cm~^ is reached within a few percent for 
N=10 layers which agrees with similar results obtained 
for the electron spectrum^. The structural, dynami- 



cal and thermodynamical properties of carbon allotropes 
were computed by Mounet et using a combination 
of density-functional theory total-energy calculations and 
density-functional perturbation theory lattice dynamics 
in the generalized gradient approximation. Very good 
agreement was found for the structural properties and 
the phonon dispersion. 

High-level ab initio calculations were carried out to 
reexamine the relative stability of bowl, cage, and ring 
isomers of C20 and (7^ by An et al.l^The total electronic 
energies of the three isomers showed a different energy 
ordering which depends on the used hybrid functional. 

Finite size planar structures which are close to the 
ground state, and in particular nanographene-like struc- 
tures were studied by Kosimov et al.^^ using energy mini- 
mization with the conjugate gradient (CG) method. The 
lowest energy, i.e. the ground state^^ configurations 
were determined for up to N=55 carbon atoms. 

Here we will investigate the normal modes for the 
ground state configuration of flat carbon clusters, also 
called nanographene as function of the number of C 
atoms in the cluster. The content of the paper is as 
follows. First (see Sec. 2), we present the main compu- 
tational approach that is used to calculate the normal 
modes. The normal modes for linear (N < 5) and ring 
(6 < < 18) clusters which are exactly one and two 
dimensional, respectively are discussed in Sec. 3. Next 
(see Sec. 4), the phonon spectrum of two dimensional 
nanographene (19< N < 55) and trigonal-shaped nan- 
odisk^^ as well as the effect of hydrogenation of the 
edge carbon atoms on the phonon spectrum are investi- 
gated. In Sec. 5, we calculate the phonon density of states 
of various clusters. By increasing the number of C atoms 
in the cluster we will show how the phonon spectrum 
of graphene appears. We also calculate the temperature 
dependence of the specific heat within the harmonic ap- 
proximation in Sec. 6. Summary (see Sec. 7) close the 
paper. 



2 



II. SIMULATION METHOD 

The Brenner second generation reactive empirical bond 
order (REBO) potential function between carbon atoms 
is used in the present work. The values for all the pa- 
rameters used in our calculation for the Brenner potential 
can be found in Ref. 31 and are therefore not listed here. 
The Brenner potential (REBO) energy is given by: 

Here Eij is the total binding energy, V^{rij) and V"^{rij) 
are a repulsive and an attractive term, respectively, 
where r^j is the distance between atoms i and j, given by 

V'ir) = f{r){l + Q/r)Ae-"^ (2) 



3 

V^{r) = f{r)J2Bne-^-'-, (3) 

n=l 

where the cut-off function f^{r) is taken from the switch- 
ing cutoff scheme 



Using molecular dynamic (MD) simulations with the 
Berendsen after corrector thermostat, we obtain the equi- 
librium configuration of stable quasi-planer clusters. We 
found that about 10^-5 x 10^ MD steps are needed to 
obtain an accuracy of 10~^-10~^ eV in the energy per 
particle. We used free boundary conditions in our simu- 
lation. 

Here, we calculate the dynamical matrix numerically 
using the finite-displacement method. We displace parti- 
cle j in the direction /3 by a small amount ±u^ and eval- 
uate the forces on every particle in the system F^- (we 

took typically u=0.003 A). Then, we compute numeri- 
cally the derivative using the central-difference formula 

dF^^, - F-^ d^E, 

- — - = lim — ■ = . (7) 



Diagonalization of the dynamical matrix yields the 
eigenvalues and the eigenvectors from which we can cal- 
culate the vibrational modes. The eigenfrequencies are 
the square root of the eigenvalues. The stability of the 
configuration can be verified because all eigenvalues have 
to be real. 
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defines the distance over which the 
function goes from one to zero and A, Q, Q^, /^n, (1 ^ 
n < 3) are parameters for the carbon-carbon pair terms. 
Here n is the type of bond (i.e. single, double and triple 
bonds). 

The empirical bond order function used in this work is 
written as a sum of terms: 



where the functions W 



and 



(5) 



depend on the local 



position and bond angles determined from their arrange- 
ment around each atom {i and j, respectively) and is 
governed by the hybridization of the orbitals around the 
atom. The first term in Eq. (5) is given by. 



[1+ E f?k{riu)G{cos{ei^u))\ 



-1/2 



(6) 



Here the angular function G{cos{Oijk)) modulates the 
contribution of each nearest neighbour and is determined 
by the cosine of the angle of the bonds between the atoms 
i-j-k. 
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FIG. 1: (Color online) Spectrum of normal modes as a func- 
tion of the number of particles in the cluster for linear N < 5 
and ring 6 < N < 18 clusters. The solid curves are the analyt- 
ical results from a simple chain (N < 5) and ring (6 < N < 18) 
model. The dashed curve indicates the breathing mode. The 
symbols are the numerical results from the present simula- 
tion where the solid (hollow) symbols are the frequencies with 
corresponding eigenvectors which are in plane (out of plane). 
The frequencies that decrease with the number of particles 
are fitted to a 1/N dependence (solid curves). 
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FIG. 2: Fitting parameter / of linear chain (N < 5) and 
ring (6 < < 18) model of particles connected by identical 
springs as a function of the number of particles. For N > 6 
the results are fitted by a simple analytical expression for even 
(full curve) and odd (dashed curve) clusters. 



III. LINEAR AND RING CLUSTERS 

The energy of the normal modes corresponding to the 
ground-state configuration for linear and ring clusters is 
shown in Fig. [l] as function of the number of C atoms in 
the clusters. The number of normal modes for N-particles 
moving in 3-dimensions is 3N. Some of these modes will 
be degenerate. There are 5 modes with uj = which 
correspond to 3 uniform translations and 2 independent 
rotations of the whole cluster, due to the translational 
and rotational invariance of the cluster. 

When clusters are exactly one dimensional, the phonon 
modes correspond to oscillations along the chain or or- 
thogonal to it. The normal modes of such a chain of N 
particles connected by springs with spring constant and 
mass m can be calculated analyticall}^ and is given by 
^ = / X "^"^ where / = 2 ^f, and 0„ = {n-l)7r/N 
{n = 1, 2, • • • , A/"). For each value of N, we determined 
/ by taking the maximum frequency obtained from our 
model equal to our numerically found result. The ob- 
tained frequencies of our linear spring chain model (solid 
curves in Fig. [T]) agree remarkably well with our numer- 
ical results. For < 5, the obtained spring constant 
is plotted in Fig. [2] and is almost constant (less than 
4% variation). As shown in Ref.^^ all the C-C bonds in 
the chain are not identical and therefore all the spring 
constants between the C-atoms are expected not to be 
identical, which is the reason why our simple model does 
not give perfect agreement with our numerical data. 

Notice that the low frequency eigenmodes are not de- 
scribed by this simple model. Some of them correspond 
to out-of-plane motion. For cluster N=2 with mode k=6 
(k counts the eigenvalues in increasing order) and for 
N=4 with mode k=12, the particles move in opposite 



direction to each other along the chain as depicted in 
Fig. [3] While the cluster N=3 with mode k=9 and the 
cluster N=5 with k=15 have alternating particles mov- 
ing in opposite direction which makes an angle with the 
chain direction (see Fig. [3|. 

The breathing mode for the linear chain is shown by 
the dotted line in Fig. [l] For N=2 the excitation cor- 
responding to the breathing mode is shown in Fig. [sj^a). 
But for N=3, this mode is similar for N=2 but where 
the direction of oscillation makes an angle with the chain 
direction and the middle particle does not move due to 
conservation of momentum. Similarly, for N=5 we found 
that the breathing mode is similar to the one for N=4 
except that the deflection makes an angle with the chain 
direction and the middle particle is at rest (see Figs.|4ja, 
b)). For mode k=6 and 7 with N=3,4,5 particles, the 
particles on the boundaries of the clusters move in the 
same direction while the remaining particles move in the 
opposite direction with different amplitudes to conserve 
total angular momentum (see Fig.[4]^c)) for both in-plane 
and out-of-plane direction. The corresponding frequency 
of these modes can be reasonably well fitted by ;^+b 
(solid curve in left panel of Fig. 1) where the fitting 
parameters are (a,b)= (1809 (±103), -234 (±27)) cm-^ 
And for mode k=8 and 9 with N=4, 5 particles, adjacent 
particles move in opposite direction with different ampli- 
tudes where the central particle is at rest in the case of 
N=5 particles (see Fig. |4jd)) for both in plane and out 
of plane direction and they can be fitted by the same 
function with (a,b)= (1783 (±103), 0.0) cm-^ 
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FIG. 3: Eigenvectors of the highest frequency modes for 
the clusters with N=2,3,4 and 5 particles for different mode 
number with frequency (a) cje = 1660 cm""*^, (b) cjg = 2092 
cm"\ (c) CJ12 2188 cm"^ and (d) cjis 2233 cm"\ 

The ground state configuration for Cn (N=6-18) con- 
sists of particles arranged in a ring. For clusters that 
are exactly two dimensional, as is the case of these ring 
structures, the phonon modes correspond to motion in 
the plane or vertical to the plane. There is no cou- 
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FIG. 4: Eigenvectors for the clusters with N=4 and 5 parti- 
cles for different mode number with frequency (a) cjio = 905 
cm""*^, (b) (jJi2 — 725 cm""*^, (c) loq — 150 cm""*^ and (d) 
(jjQ = 466 cm~-^. 
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pling between the in-plane and the out-of-plane motion. 
We can construct a simple model of equidistant par- 
ticles arranged in a ring, where neighbor particles are 
connected by a spring. This model can be viewed as 
the linear chain model with periodic boundaries. For 
fixed radius of the ring this results in the following 
eigenfrequencieP^ = / x sin{k/2)^ where k = 27rn/N^ 
(n = 1, 2, • • • , A^). As before we determined / such that 
the maximum frequency obtained from the model equals 
the numerical result. This fitted value is shown in Fig. |2] 
which can be described approximately by the function 
f{N) = ax where: a=3200 cm~-^ and i) for even 

N clusters b = -0.295 (± 0.001), and c = -0.395 (± 0.002); 
and ii) for odd N clusters b = -0.276 (± 0.001), and c = 
-0.369 (± 0.001). Several of the numerical frequencies are 
reasonably well described by this simple model. Notice 
that the eigenfrequencies exhibit a clear discontinuous 
behavior when a cluster changes from a linear into a ring 
configuration. As an example we show in Fig. [5] some of 
the interesting eigenmodes of a cluster with N=18 parti- 
cles which are arranged in a circle and is therefore pure 
two dimensional. The mode k=23 is the breathing mode 
while mode k=38 corresponds to radial out of phase ra- 
dial oscillation of nearest neighbor particles. The modes 
k=53, and 54 correspond to pure angular oscillations of 
the particles. 

It is interesting to note that for an even number of par- 
ticles, there is a nearly common frequency (i.e. N=8 (545 
cm-i), N=10 (543 cm-^), N=12 (546 cm-^), N=14 (547 
cm-i), N=16 (549 cm-^), N=18 (550 cm-^)) which is 
almost independent of N (see Fig. fl]). This mode corre- 
sponds to out of phase oscillations of nearest-neighbor 
atoms (see Fig. |5|b) the k=38 mode for N=18 and 
Fig. [6]^ a) the k=24 mode for N=12). But for odd ring 
clusters i.e. N=13, the corresponding normal mode corre- 



FIG. 5: Eigenvectors for a cluster with N=18 particles for 
different mode number with frequency (a) (jJ23 = 398 cm~^, 
(b) UJ38 550 cm"\ (c) UJ53 2224 cm"^ and (d) 2259 
cm""*^. 



sponds to nearest-neighbor atoms oscillating out of phase 
with different amplitudes except for the two neighbor 
atoms that oscillate in phase (see Fig. |6jb) k=23 (471 
cm~^)) which results in a lower frequency by about 20 % 
as compared to the corresponding even N ring clusters. 

A breathing mode exists for all the ring clusters as 
shown in Fig. [6) The frequencies for N=8, 12, 13, 16, 
and 18 are 837 cm~^, 585 cm~^, 543 cm~^, 447 cm~^, 
and 398.25 cm~^, respectively, and have a clear depen- 
dence on N. The simple spring model^^ predicts that this 
breathing mode has the frequency uj={7r x f)/N which 
remarkably agrees with our numerical results. 

The eigenmodes corresponding to uJrnax have similar 
vibrational modes for even number of particles. We plot- 
ted the normal modes of ujmax of cyclic structures for 
N=12 {(jJsQ = 2183 cm~^) as shown in Fig. [6|e) and they 
correspond to dipole-type of oscillations between nearest 
neighbor particles while for odd number of ring struc- 
tures i.e. N=13 (cjag = 2187 cm~^ in Fig. [6|f)), similar 
dipole-type of oscillations between nearest neighbor par- 
ticles are found but with decreasing magnitude towards 
opposite sides. 

Now we investigate the eigenmodes for out-of-plane vi- 
brations for ring clusters. For N=10 and mode k=10, the 
normal mode corresponds to a bending mode while for 
mode k=16, the normal mode corresponds to a sinusoidal 
type of motion. The higher modes k=17 and k=19 have 
nearest-neighbor atoms oscillating in opposite direction 
perpendicular to the ring-plane with different wavelength 
along the ring as seen in Fig. [7[ For even ring structures 
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FIG. 6: Normal modes for the out-of-phase (top panels), 
breathing (middle panels) and Umax (bottom panels) for the 
clusters N=12 (left panel) and N=13 (right panel). 



TABLE I: Fitting parameters for the eigenfrequencies of out- 
of-plane motion shown by the solid curves in Fig.l for C(;<550 




FIG. 7: Eigenvectors for out-of-plane motion for a cluster 
with N=10 particles for different mode number with frequency 
(a) cjio = 165 cm~\ (b) cjie = 482 cm~\ (c) cjit = 482 cm~^ 
and (d) c^ig — 536 cm~^. 
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FIG. 8: 1/Dc-c (dashed curve) and uJraax (solid curve for 
even and odd N) as a function of the number of particles. 



2782 (±48) -107 (±5) 
4258 (±120) -97 (±11) 
5405 (±251) -62 (±20) 
5686 (±59) 25 (±5) 
5814 (±499) 102 (±33) 
4942 (±589) 220 (±36) 



the highest mode for out-of-plane vibration is similar to 
the one shown in Fig. [7|^d). The out-of-plane vibration 
can be reasonably well fitted by ± 6 (solid curves in 
Fig. 1 for a;<550 cm~^) where the fitting parameters a 
and b are listed in Table [ij The eigenmodes belonging to 
each solid curve correspond to the same type of mode. 



From Fig.[T] it is interesting to note that the maximum 
frequency in the excitation spectrum for ring clusters, 
on the average, slowly increases with N. This can be ex- 
plained by calculating the phonon spectrum of an infinite 
system. In Fig. [8] we show that the minimal interparticle 
distance decreases slowly with the number of C-atoms 
in the cluster. As a consequence, the maximum value 
of the wave vector /c ~ tt/Zq (^o is the mean distance be- 
tween the particles) and also the corresponding frequency 
will increase weakly with cluster size. The inverse of the 
inter-carbon distance could be fitted to (dashed curve in 



Dc-G Dol + bN 



6 



with the fitting parameters: -^=7.9 nm ^, a = -0.316 
(± 0.002), and b = -0.331 (± 0°.002). We are able to fit 
a curve through the maximum frequency of these ring 
clusters using (solid curves in Fig. [8| 



1 + a7V 



(9) 



where u;o=3200 cm ^, and: i) for even N clusters a = 
-0.295 (± 0.001), and b = -0.395 (± 0.002); and ii) for 
odd N clusters a = -0.276 (± 0.001), and b = -0.369 (± 
0.001). 



IV. NANOGRAPHENE {N > 18) 

The clusters with > 18 have an inner structure and 
their configurations have been investigated in detail in 
Ref. |27| Their ground state configuration can be clas- 
sified in three groups: 1) nanographene clusters consist- 
ing of only hexagons, 2) clusters with pentagon on the 
boundary, and 3) bowl shaped configurations that have 
typically pentagons in the inner part of the cluster. For 
the latter one < z'^ >^ 0, where z is the position co- 
ordinate of the cluster along the out-of-plane direction 
and they are found for N=20, 28, 38 and 44 which are 
buckled-like structures. In these clusters one pentagon is 
surrounded by five hexagons. The normal mode oscilla- 
tions in-plane and out-of-plane are now coupled, i.e. the 
normal mode oscillations are 3-dimensional and some of 
the interesting ones are shown in Fig. [9] for N=20 and 
mode k=7, 9, 11 and 60. For mode k=7, the opposite 
atoms on the boundary of the cluster vibrate in the same 
direction whereas for k=9, the normal mode corresponds 
to a bending mode of the cluster. But, for the higher 
mode k=ll, the corner atoms show mixed type of os- 
cillations while for large value of k, the normal modes 
correspond to in-plane oscillations of the C-atoms. 

The normal mode frequencies for A/" > 19 are plotted in 
Fig. 10 For N=22 and 39, a heptagon is on the boundary 
surrounded with pentagon and hexagon. In such a case 
local modes are found where only particles close to the 
defects (pentagon and heptagon) oscillate. If there is 
one pentagon on the boundary i.e. N=49, and 51, the 
local modes with larger amplitudes are on the opposite 
side of the clusters while for N=53 where a pentagon is 
surrounded by four hexagons, the local modes with larger 
amplitudes are found only near the defect. For clusters 
which are symmetric and without defect i.e. N=54, we 
found that the local modes are situated on the boundary 



as shown in Fig. 11 



A special subgroup of clusters consists of hexagonal 
and trigonal shaped nanographene. The clusters with 
size N=24 and N=54 are planar graphene structures 
which are hexagonal shaped with zigzag edges. These 
clusters have a close-packed structure, which consist 
purely of hexagons. Let us consider the cluster with 
N=54 (see Fig. 12), the mode k=34 corresponds to the 
breathing mode. Mode k=58 exhibits circular motion of 
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FIG. 9: Eigenvectors for a cluster with N=20 particles for 
different mode number with frequency (a) ur = 124 cm~^, 
(b) cjg — 161 cm""*^, (c) cjii — 233 cm""*^ and (d) cjeo — 1772 
cm~^. 
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FIG. 10: (Color online) Normal mode spectrum as a function 
of the number of C atoms for 19 < N < 55. The modes with 
eigenvectors in-plane are shown in black, those with out-of- 
plane eigenvectors in red and mixed type eigenvectors in blue. 



atoms arranged in hexagons near the 6 corners of the 
hexagonal disk. The next mode k=59 is similar but now 
hexagons are displaced by an angle of 30° and are cen- 
tered in the middle of the sides of the hexagonal disk. 
Notice that for k=58 the 6 rotational oscillations are in 
the same direction while for k=59 they alternate, i.e. 
vortex/anti- vortex like arrangements. For higher modes 
k=139, only the inner particles participate in the normal 
mode oscillations. For k=148, a dipole type of oscilla- 
tions is found of nearest neighbor C-atoms arranged in 
a shell around the middle between the center and the 
perimeter. For higher frequency only the outer particles 
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FIG. 12: Eigenvectors for a cluster with N=54 particles for 
different mode number with frequency (a) (JJ34 = 313 cm~^, 
(b) cjss = 531 cm~-^, (c) cjsg = 534 cm~-^, (d) cjro = 593 cm~-^, 
(e) CJ139 = 1589 cm"\ (f) cji48 = 1624 cm"\ (g) CJ149 = 1628 
cm~^, (h) (jji^o — 1638 cm~^ and (i) cji62 = 1914 cm~^. 



FIG. 11: Eigenvectors for highest mode number for clusters 
with N=22, 39, 49,51, 53, 54 particles with frequency (a) 
CJ66 = 1944 cm"\ (b) cjut = 1926 cm"\ (c) CJ147 = 1907 
cm~-^, (d) cjiss = 1903 cm~-^, (e) cjisg = 1952 cm~-^ and (f) 
UJIQ2 — 1914 cm~^. 



oscillate while the inner particles exhibit only very small 
displacements as shown in Fig. [l2j 

From Fig. [To] we notice that there is a region along the 
frequency axis with a higher density of normal modes. 
Examples of the displacements of such modes are shown 
in Fig. [13] for N=54. Mode k=27 corresponds to the ex- 
citation of a vortex/antivortex pair. Mode k=33 consists 
of rotational oscillations around the hexagonal corners 
of the nanodisk, while for k=35 these rotational motions 
are centered around the middle of the hexagonal sides. 
Notice that mode k=45 consists of a central large vor- 
tex motion surrounded by 6 anti- vortex type of motions 
situated closer to the edge of the nanodisk. 

The nanodisk cluster with N=46 carbon atoms has 
a metastable configuration consisting of closed zigzag 
edges and a trigonal-shaped structur^^. As shown in 
Fig. [14] the mode k=18 consists of an asymmetric vor- 
tex/antivortex pair while mode k=29 is a clear breathing 
mode. The modes k=24 and k=32 show three rotations 



situated close to the three corners of the trigonal nan- 
odisk, but notice that the rotation direction and the po- 
sition of the center of the rotation is not always the same 
for both modes . 

We plotted the average distance between the C-C 
atoms as a function of the number of particles N. The 
average radius increases linearly with the chain length for 
(N=3-5) which could be fitted to Dc-c = Dq ^ a x N 
where Dq = 0.13 nm and a=0.00055 (± 0.00009) nm (see 
solid curve in Fig. 15). For ring clusters it decreases ex- 
ponentially as Dc-c = Dq -\- a X exp{b x N) (see solid 



curve in Fig. [15| where the parameter I^o = 0.134 nm, 
a = 0.039(±0.001) nm and b = -0.337(±0.007) as the 
number of particles increases. For N > IS the average C- 
C distance fluctuates as function of N around the average 
value 0.141 nm which compares with 0.142 nm for the C- 
C distance in bulk graphene. The clusters with pentagon 
and Stone- Wales defects inside the clusters show larger 
average C-C distance than the pure hexagonal structures 
while the clusters with pentagons on the boundary (i.e. 
N=26, 31, 34, 36, 41, 46, 51 and 53) show the highest 
average distance as shown in Fig. [T5| 



For completeness we also plotted (in Fig. 16) the en- 
ergy spectrum of the normal modes for ground-state H- 
passivated Cat (N= 11-55) cluster^!. Notice that dif- 
ferent from Fig. [To] (i.e. non H-passivated graphene) 
there is a region with enhanced density of modes, i.e. 
around 2900 cm~^. The corresponding modes are con- 
nected to oscillations of the C-H atoms. This is illus- 
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FIG. 13: Eigenvectors for the cluster with N=54 particles 
for four different values of the mode number k with frequency 
(a) 0027 244 cm"\ (b) CJ33 302 cm"\ (c) UJ35 314 cm"^ 
and (d) CJ45 = 408 cm~X 



FIG. 14: Eigenvectors for a cluster with N=46 particles for 
different mode number with frequency (a) c^is = 185 cm~^, 
(b) CJ24 244 cm"\ (c) UJ29 288 cm"^ and (d) CJ32 322 
cm""*^. 



trated in Fig. [T7| for H-passivated nanographene clusters 
with N=53 and 54. In the optical region which is a re- 
gion along the frequency axis with higher density of nor- 
mal modes, for N=53 with k=172 the corresponding local 
normal modes with larger amplitudes are on the opposite 
side of the defect while for N=54 with k=187, only the 
inner C-atoms oscillate with larger amplitude (Figs. 17 
(a, b)). The highest frequency mode of the N=53 and 
54 clusters (Figs 



17 



(c, d)) is a local mode near the lo- 
cal defect where mainly outer C-H atoms of the armchair 
C-atoms oscillate. 



V. PHONON DENSITY OF STATES 

In order to compare the normal modes with those 
of graphene we calculate the phonon density of states 
(PDOS). Furthermore, we calculated the PDOS of ex- 
actly trigonal and hexagonal two dimensional and de- 
fective clusters and analyze what is the effect of defects 
on the density of states. Because of the discreteness of 
the normal modes frequency spectrum we introduced a 
Gaussian broadening 



3N 



(10) 



where the summation is over all normal modes, uji is the 
normal mode frequency of the i^^ mode, broadening is 
chosen Suj = 30 cm~^ and N is the total number of atoms 
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FIG. 15: The average distance between C-C atoms as a 
function of the number of particles in the graphene cluster. 
The error bars indicate the range of C-C distances within each 
cluster. 



in the cluster. We notice that the introduction of a pen- 
tagon in the N=53 cluster as compared to the perfect 
hexagon lattice in N=54 introduces more high frequency 



modes as shown in Fig. 18 which are due to modes local- 
ized around the defect. 



Fig. 19 shows the phonon density of states (PDOS) of 
armchair hexagonal and zigzag trigonal clusters. Notice 
that several of the peaks in the PDOS of the hexagonal 
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FIG. 16: (Color online) Normal mode spectrum as a function 
of the number of C particles for 11 < N < 55 H-passivated 
clusters. The modes with eigenvectors in plane are shown in 
black, those having out-of plane eigenvectors in red. 
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FIG. 17: Eigenvectors for a cluster with N=53, and 54 with 
H-passivation for k=172 for N=53 and k=187 for N=54, and 
for ujmax corresponding to frequency (a) = 1693 cm~^, 
(b) cjis? = 1698 cm~-^, (c) 0)210 = 2925 cm~-^ and (d) 00222 = 
2915 cm-\ 



clusters are split into two peaks in the PDOS of the trig- 
onal shaped cluster. This is a consequence of the reduced 
rotational symmetry of the trigonal cluster. 



We compare in Fig. [20] the phonon density of states of 
clusters N=55 and 1600 with those for graphene^^ which 
shows nicely the convergence of the phonon spectrum of 
nanographene to bulk graphene. For N=1600 we almost 
recover the theoretical results of graphene^^ . Notice that 
for small clusters there is a pronounced PDOS for a; = 
which is due to the relative importance of the uj = 
translational and rotational motion. 
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FIG. 18: (Color online) The phonon density of states for 
clusters with N= 53 and 54 carbon atoms. 
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FIG. 19: (Color online) The phonon density of states for 
clusters with N=42 (hexagonal with armchair) and 46 (trigo- 
nal) carbon atoms. 



The effect of H-passivation of the edge atoms of the 
cluster is shown in Fig. [2l] for N=53 and 55. The C- 
H bonds introduce high frequency modes in the N=53 
and N=55 clusters as compared to the non H-passivated 
N=53 and N=55 clusters. The H-passivated clusters ex- 
hibit pronounced peaks around 1700 and 2900 cm~^ 
which correspond to the optical region and the C-H atom 
frequency region, respectively. We also notice that the 
introduction of a pentagon in the N=53 cluster as com- 
pared to the perfect hexagon lattice for N=55 introduces 

similar to 



more high frequency modes as shown in Fig. [21 



18) which 



the case for non H-passivated clusters (see Fig. 
are due to modes localized around the defect. 

The modes can have a shear-like or a compression-like 
character. The compressional and shear properties can 
be extracted from the divergence and rotation of the ve- 
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FIG. 20: (Color online) The phonon density of states for 
clusters with N=55, 1600 and perfect graphene (bottom panel 
is taken from Ref . ,35j . 
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FIG. 21: (Color online) The phonon density of states for 
H-passivated clusters with N=53 and 55. 



locity field. Here, we will associate a single number for 
the shear-like and compression-like character of the dif- 
ferent modes by calculating the squared average of the 
divergence V.v and the vorticity (V xv)^ of the velocity 
field, following the approach of Ref. l36l 

The z component of the rotation ^^(/c) = ez.rot^ik) 
and the divergence ^d{^) = div^(k) of the field of the 
eigenvectors of mode k are 



1 ^ 
N ^ 



iik), 



(11) 



1 ^ 



(12) 



where the values of "^d^ik) and ^r,i(^) for the particle 
are, respectively, given by 

1 ^ 

"^dAk) = -J - ^j)'[Mk) - - r,f , (13) 



1 ^ 

-^(rl-r,Ox[l,(fc)-A,(fc)]/|r-;-r,f. (14) 



J 



Here, j and J denote the index and the number of neigh- 
boring particles of particle i, respectively, fj is the posi- 
tional coordinate of a neighboring particle and Ai{k) is 
the eigenvector of particle i for mode k. Note also that 
we calculate the squared average over all the particles 
because the simple spatial average is of course zero. 

We plot ^d(^) and ^r(^) as a function of the mode 
number k for cluster N=39 which contains a 5-7 defect, 
cluster N=49 with pentagon defect on the boundary, the 
pure symmetric hexagonal cluster with N=54 and for a 
large cluster having N=398 particles. In general, the 
lower eigenfrequency spectrum corresponds to rotational 
type of excitations which are vortex- ant i vortex like ex- 
citations for large clusters (see Fig. 13). Vortex excita- 
tion is only expected for sufficiently large clusters, be- 
cause the velocity of dissipation of the vortex energy is 
inversely proportional to R^, where R is the radius, which 
increases with increasing cluster size. In the second half 
of the spectrum the divergence "^dik), which corresponds 
to compression-like modes, is appreciably different from 
zero and we have mixed modes that have both a shear-like 
and a compression-like component. For symmetric clus- 
ters and defective clusters with pentagon on the bound- 
ary, the highest modes consist of only shear-like modes 
(see Figs. 11 ^c, d, f)). For cluster N=39 with a 5-7 defect, 
compression-like modes appear for the highest frequency 
modes (see Fig. [iTJb)). 

The heat content of a system is directly related to the 
specific heat whose lattice contribution is determined by 
the phonons i.e. the eigenmodes of the system, while 
electronic contribution to it can be neglected even at a 
few Kelvins^^. The specific heat Cp per unit mass, at 
temperature T is within the harmonic approximatioiP^ 
given by the following expression: 



J 3N 

ku sr^ 
MN ^ 

k=l 



2kBT 



csch 



2kBT 



(15) 



where the summation is over all normal modes, kp is 
the Boltzmann constant, M is the mass of each carbon 
atom and N is the total number of atoms in the cluster. 
From Eq. (15), the specific heat depends sensitively on 
the characteristics of the phonon spectrum and on its 
vibrational density of states. 
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FIG. 22: (Color online) Distribution of the rotation and 
the divergence as function of the mode number k for N=39, 
49, 54 and 398 particles. 



TABLE II: Fitting parameters for the specific heat (Eq. (16)) 
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FIG. 23: (Color online) Phonon contribution to the constant 
pressure heat capacity of clusters with N=5, 18, 54 and 398 at 
low temperature. The curves are linear fits to the numerical 
results. In the inset we show the N-dependence of the slope of 
the low temperature linear T-dependence of the specific heat. 
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The numerical results for N=5, 18, 54 and 398 are 
It may be noted that at low tem- 
=5, 18, 54 and 398 have a 



shown in Fig. 23 
perature, cluster sizes with N= 
specific heat Cp which exhibits a clear linear dependence 
in T. The slope of the low temperature Cp behavior is 
shown in the inset of Fig. [23] as a function of N which 
can be fitted to: a ± b x exp(c x N) (solid curve) where 
a=3.77 J kg-i K-\ b=3.87 (±1.23) J kg-^ K'^ and 
c=-0.05 (±0.01) K-^ 



The high temperature behavior is shown in Fig. 24 For 
large temperatures all results approach the value 2078 
J kg~^ which is in good agreement with the bulk 
graphene value given in Refs) ^^*^^ *. 

We fitted the numerical results using the formula 



Cp = a ± 6 X exp{c x T), 



(16) 



where a=2078.0 J kg~^ and b and c are fitting pa- 
rameters. After fitting the curve, we obtained the pa- 
rameters listed in Table [ll} Notice that \b\ increases with 
N while IcI decreases. 
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FIG. 24: (Color online) Phonon contribution to the constant 
pressure heat capacity of clusters with N=5, 18, 54 and 398 at 
high temperature. The numerical results approach the fitted 
value 2078 Jkg-^K'K 



VI. SUMMARY 

Atomistic simulations were performed using the Bren- 
ner second-generation reactive bond order (REBO) in- 
teratomic potential function in order to study the vibra- 
tional properties of nanographene. In the present work 
we investigated the frequency of the normal modes as 
a function of the number of particles for nanographene 
Cat (2 < < 55) and H-passivated carbon clusters 
(11 < N < 55). We also presented a simple model for the 
eigenfrequencies for linear chain and ring structures and 
obtained analytical results that are in good agreement 
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with our numerical results. The phonon density of states 
of different clusters are compared with theoretical results 
for graphene and found to be in good agreement for the 
acoustic and optical phonon density in the case of large 
clusters. Using the eigenfrequencies of the normal modes 
we calculated the specific heat within the harmonic ap- 
proximation for different size carbon clusters and found 
that in the high temperature limit they approach the bulk 
value as given in Refs.^^'^. 
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